Imports

library(evd); 
library(evdbayes);
library(coda);
library(ismev);
Lade nötiges Paket: mgcv
Lade nötiges Paket: nlme
This is mgcv 1.8-24. For overview type 'help("mgcv-package")'.
library(xts);
Lade nötiges Paket: zoo

Attache Paket: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

Loading the Data

load("../data/CAPE_Minder_Rychener_Malsot.RData")
load("../data/NINO34.RData")
load("../data/SRH_Minder_Rychener_Malsot.RData")

Generate PROD

prod = sqrt(cape)*srh

** Create Time Series Objects **

dates <- seq.Date(as.Date("1979-1-1"), as.Date("2015-12-31"), by=1)
feb29ix <- format(as.Date(dates), "%m-%d") == "02-29"
dates <- dates[!feb29ix]
prod_ts = xts(prod, order.by = rep(dates, each=8))

Beginning the Analysis

month_names = c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
get_monthly = function(x) {
  output = list()
  len = nrow(x)
  
  # Get month of element
  month = time(x)
  month = gsub("....-", "", month)
  month = gsub("-..", "", month)
  monthlist = unique(month)
  for (i in 1:12) {
    output[[i]] = x[month == monthlist[i],]
  }
  names(output) = month_names
  return(output)
}
monthly_max = get_monthly(apply.monthly(prod_ts, max))
r = 2
r_monthly_max = get_monthly(apply.monthly(prod_ts, function(x) order(x, decreasing=T)[1:r]))
# get the monthly maxima of prod
m1 = as.data.frame(apply.monthly(prod_ts, max))$V1;
# produce the gumbel qq plot
gumbel_qq = function(x, title="") qqplot(qgumbel(c(1:length(x))/(length(x)+1)),
                                         x,
                                         main = paste("Gumbell Q-Q Plot", title),
                                         xlab = "Theoretical Quantiles", 
                                         ylab = "Sample Quantiles") ;
gumbel_qq(m1)

#qqplot(qgumbel(c(1:length(m1))/(length(m1)+1)),m1,main = "Gumbell Q-Q Plot",xlab = "Theoretical Quantiles", ylab = "Sample Quantiles") ;

The qq plot doesn’t fit very well, especially in the lower tail. This is likely due to seasonal dependence.

Fitting GEV to the entire data

# fit gevd with MLE and produce diagnostic plots
fitmax.MLE<-fgev(m1,method="Nelder-Mead")
par(mfrow=c(2,2))
fitmax.MLE

Call: fgev(x = m1, method = "Nelder-Mead") 
Deviance: 9272.599 

Estimates
      loc      scale      shape  
7.519e+03  7.013e+03  4.757e-03  

Standard Errors
      loc      scale      shape  
421.90201  331.43749    0.05347  

Optimization Information
  Convergence: successful 
  Function Evaluations: 171 
plot(fitmax.MLE)

Poor fit, probably because the distribution isn’t stationary. This is especially visible in the Probability plot, in which the confidence band is surpassed, indicating a poor fit.

# fit gevd with Bayesian Techniques
# use the previous outputs (rounded) as initial values (use a different shape)
init<-c(1.6e4,4e3,0.1)
# arbitrary priors
mat <- diag(c(10000,10000,100)) 
pn <- prior.norm(mean=c(0,0,0),cov=mat)
# proposal standard deviation (found by trial and error to get 20%<acceptance rate<40%)
psd<-c(500,0.03,0.02)
# draw 3k samples from markov chain
nit = 30000
thinning = 300
post<-posterior(nit, init=init,prior=pn,lh="gev",data=m1,psd=psd)
# diagnostic plots
MCMC<-mcmc(post[1:nit %% thinning == 0, ], thin=300) 
plot(MCMC) 

attr(mcmc(post),"ar")
            mu sigma   xi total
acc.rates 0.24  0.71 0.70  0.55
ext.rates 0.00  0.00 0.02  0.01
#MCMC_stab <- mcmc(post, thin=50, start=1000)
acf(mcmc(post[1:nit %% thinning == 0, ]))

We observe that there seem to be no substantial issues from the autocorrelation plots.

apply(mcmc(post[1:nit %% thinning == 0, ]),2,mean)
           mu         sigma            xi 
  214.8136408 11344.8681059    -0.1940404 
apply(mcmc(post[1:nit %% thinning == 0, ]),2,sd)
          mu        sigma           xi 
100.50420904 741.34912578   0.03950923 

** Fit with MLE for months separately**

#monthly_fits = lapply(monthly_max, 
#                      function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
monthly_fits = list()
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
  print(i)
  
  if (i %in% error_cases) {
    monthly_fits[[i]] = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             std.err = FALSE)
  }
  else {
    monthly_fits[[i]] = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead")
  }
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
optimization may not have succeeded
[1] 10
[1] 11
[1] 12
names(monthly_fits) = names(monthly_max)
gumbel_qq(data.frame(monthly_max[[9]])[,1], "September")

gumbel_qq(data.frame(monthly_max[[12]])[,1], "December")

get_se = function(x, ix) {
  if (is.null(x$std.err)) 0
  else x$std.err[ix]
}
mle_loc = unlist(lapply(monthly_fits, function(x) x$estimate[1]))
mle_loc_se = unlist(lapply(monthly_fits, get_se, 1))
mle_scale = unlist(lapply(monthly_fits, function(x) x$estimate[2]))
mle_scale_se = unlist(lapply(monthly_fits, get_se, 2))
mle_shape = unlist(lapply(monthly_fits, function(x) x$estimate[3]))
mle_shape_se = unlist(lapply(monthly_fits, get_se, 3))
plot_w_err = function(x, y, se, title = NULL) {
  max_ix = which.max(y)
  min_ix = which.min(y)
  plot(x, y,
       ylim = c(y[min_ix] - se[min_ix], y[max_ix] + se[max_ix]),
       main = title)
  arrows(x,y-se,x,y+se, code=3, length=0.02, angle = 90)
}
plot_w_err(1:12, mle_loc, mle_loc_se, "Location Parameter as Estimated with Likelihood")

plot_w_err(1:12, mle_scale, mle_scale_se, "Scale Parameter as Estimated with Likelihood")

plot_w_err(1:12, mle_shape, mle_shape_se, "Shape Parameter as Estimated with Likelihood")

** Fit with Bayesian Methods for months separately**

# Fits GEV distribution with bayesian method for given parameters
bayes_fitter = function(x, 
                        init = c(1e3, 1e3, 0.1), # Initial values
                        mat = diag(c(10000,10000,100)),
                        psd = c(500,0.1,0.1), # Proposed SDev
                        nit = 3000, # Nb Iterations
                        thinning = 50, # Thinning Factor
                        do_diagn = FALSE, # Bool whether to show diagnostic plots
                        do_autoreg = FALSE, # Bool whether to show autoreg plots
                        seed = 42 # Seed 
                        ) {
  set.seed(seed)                
  pn = prior.norm(mean=c(0,0,0),cov=mat)
  post = posterior(nit, init=init, prior=pn, lh="gev", data=x, psd=psd)
  
  if(do_diagn) {
    MCMC<-mcmc(post[1:nit %% thinning == 0, ]) 
    plot(MCMC) 
  }
  if(do_autoreg) {
    acf(mcmc(post[1:nit %% thinning == 0, ]))
  }
  list(posterior = post, 
       acceptance_rate = attr(mcmc(post),"ar"))
}
# Iteratively fits GEV with bayesian methods, until the fit has 
# acceptable acceptance rates (i.e. 0.2 < AR < 0.4). If the AR is too high, 
# the proposed SDev for the parameter is multiplied with 1.5. If it's too small,
# the proposed SDev is divided by 2. This is repeated until the acceptance rate
# is good for all parameters, or max_it is reached. Then, a final model is fitted
# with more iterations. 
bayes_fitter_param_search = function(x,
                                     psd_init = c(500,0.1,0.1), # Initial proposed SDev
                                     nit_full = 3000, # Nb Iterations for final model
                                     nit_search = 150, # Nb Iterations for param search
                                     do_diagn = FALSE, # Bool whether to show diagnostic plots
                                     do_autoreg = FALSE, # Bool whether to show autoreg plots
                                     max_it = 20,
                                     ... # Additional params passed to bayes_fitter
                                     ) 
{
  # Iterate until desired acceptance rate is obtained
  cont = TRUE
  psd = psd_init
  it = 0
  while(cont) {
    it = it+1
    if (it > max_it) {
      warning("The ")
    }
    fit = bayes_fitter(x, psd=psd, nit=nit_search, do_diagn=FALSE, 
                       do_autoreg=FALSE,...)
    acc_rates = fit$acceptance_rate[1, 1:3]
    
    too_high = acc_rates > .4
    too_low = acc_rates < .2
    
    if (all(!too_high) && all(!too_low)) { # All acceptance rates lie within threshold
      cont=FALSE
    } else if (it > max_it) { # max_it is reached
      cont=FALSE
      warning("max_it was reached")
    } else { # Correct values which have wrong threshold
      psd[too_high] = psd[too_high] * 1.5
      psd[too_low] = psd[too_low] / 2
    }
  }
  
  # Fit final model
  bayes_fitter(x, psd=psd, nit=nit_full, do_diagn=do_diagn, 
               do_autoreg=do_autoreg, ...)
  
}
                                     
monthly_bayes_fit = lapply(monthly_max, bayes_fitter_param_search, do_diagn = TRUE, 
                           do_autoreg = FALSE,
                           psd = c(500,0.3,0.3), nit_full=30000, nit_search = 30000,
                           thinning = 300)

TODO-> R largest fit

PART 2 First, we check if the location parameter depends on time using a likelyhood ratio test

#monthly_fits = lapply(monthly_max, 
#                      function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
ratios = list()
trend = 1:length(as.data.frame(monthly_max[[12]])$V1)
trend = (trend-mean(trend))/sd(trend) # scale and center covariates as recommended
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
  print(i)
  
  fit_const = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             std.err = FALSE)
  fit_dependant = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             nsloc = trend,
                             std.err = FALSE)
  
  ratios[[i]] = fit_const$dev-fit_dependant$dev 
}
[1] 1
optimization may not have succeeded
[1] 2
optimization may not have succeeded
[1] 3
[1] 4
optimization may not have succeeded
[1] 5
optimization may not have succeeded
[1] 6
[1] 7
[1] 8
[1] 9
optimization may not have succeeded
[1] 10
[1] 11
optimization may not have succeeded
[1] 12
names(ratios) = names(monthly_max)
chi_95level = qchisq(1-0.05/12,1)
plot(unlist(ratios),main="95% confidence test for time independance, Bonferroni multiple Testing", xlab="Month",ylab="Likelyhood Ratio Statistic")
abline(a=chi_95level,b=0,col="red")

Now, let’s check for independance from ENSO

#monthly_fits = lapply(monthly_max, 
#                      function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
ratios = list()
# split nino data into months
n = nino34
dim(n)=c(12,length(as.data.frame(monthly_max[[12]])$V1))
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
  print(i)
  trend = n[i,]
  trend = (trend-mean(trend))/sd(trend) # scale and center covariates as recommended
  fit_const = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             std.err = FALSE)
  fit_dependant = fgev(as.data.frame(monthly_max[[i]])$V1, 
                             method = "Nelder-Mead",
                             nsloc = trend,
                             std.err = FALSE)
  
  ratios[[i]] = fit_const$dev-fit_dependant$dev 
}
[1] 1
optimization may not have succeeded
[1] 2
optimization may not have succeeded
[1] 3
[1] 4
optimization may not have succeeded
[1] 5
[1] 6
optimization may not have succeeded
[1] 7
[1] 8
optimization may not have succeeded
[1] 9
optimization may not have succeeded
[1] 10
[1] 11
[1] 12
names(ratios) = names(monthly_max)
chi_95level = qchisq(1-0.05/12,1)
plot(unlist(ratios),main="95% confidence test for independance from ENSO, Bonferroni multiple testing", xlab="Month",ylab="Likelyhood Ratio Statistic")
abline(a=chi_95level,b=0,col="red")

Another method is the chi plot:

# plot the chi plot for dependance analysis
trend = 1:length(as.data.frame(monthly_max[[12]])$V1)
n = nino34
dim(n)=c(12,length(as.data.frame(monthly_max[[12]])$V1))
for (i in 1:length(monthly_max)) {
  nino = n[i,]
  m_data=as.data.frame(monthly_max[[i]])$V1
  dat.m1_month = cbind(m_data,trend);
  dat.m1_nino = cbind(m_data,nino);
  par(mfrow=c(2,2))
  chiplot(dat.m1_month,main1 = "Chi Plot Time",main2 = "Chi Bar Plot Time");
  chiplot(dat.m1_nino,main1 = "Chi Plot ENSO",main2 = "Chi Bar Plot ENSO");
}

PART 3 We will now analyse temporal clustering of extremes. For this, we will use the exiplot function from the evd library.

# define a function for getting the extremal indices for each month for a given threshold
monthly_eindex <- function(data, threshold_p, r=0){
  ei = list()
  for (i in 1:length(data)) {
    threshold = quantile(as.data.frame(data[[i]])$V1, threshold_p)
    ei[[i]]=exi(as.data.frame(data[[i]])$V1, threshold, r)
  }
  names(ei) = names(data)
  return(ei)
}
ei = monthly_eindex(get_monthly(prod_ts), 0.95)
plot(unlist(ei), main="Extremal Index by Month, 95%-Quantile as Threshold", xlab="Month", ylab="Extremal index")

We observe that the extremal index is 0.25-0.45, we can therefore conclude that we have strong dependance of extremes, with the limiting mean cluster size being roughly from 2 to 4. The clustering has no effect for estimaters based on the (monthly) maximum, but the r largest estimator is influenced by it.

PART 4 First, let’s estimate the 10 year return level using point process approach

monthly_fits_pp = list()
monthly_data = get_monthly(prod_ts)
#error_cases = c(1,2,3,4,5,6,7,8,9,10,11,12)
month_days = c(31,28,31,30,31,30,31,31,30,31,30,31)
for (i in 1:length(monthly_max)) {
  print(i)
  threshold = quantile(as.data.frame(monthly_data[[i]])$V1, 0.95)
  
  if (i %in% error_cases) {
    monthly_fits_pp[[i]] = fpot(as.data.frame(monthly_data[[i]])$V1,
                             threshold = threshold,
                             model="pp",
                             npp = month_days[i]*8,
                             cmax = TRUE,
                             r = 1,
                             std.err = FALSE,
                             method = "Nelder-Mead")
  }
  else {
    monthly_fits_pp[[i]] = fpot(as.data.frame(monthly_data[[i]])$V1,
                             threshold = threshold,
                             model="pp",
                             npp = month_days[i]*8,
                             cmax = TRUE,
                             r = 1,
                             method = "Nelder-Mead")
  }
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
names(monthly_fits_pp) = names(monthly_data)
for(i in 1:12){
  par(mfrow=c(2,2)) 
  plot(monthly_fits_pp[[i]])
}

The fit in february has completely failed, and the others are not very good either

We will still estimate the return levels:

return_level = function(x,period=20){
  p = 1/period
  loc = x$estimate[[1]]
  scale = x$estimate[[2]]
  shape = x$estimate[[3]]
  level = loc + scale*(((-log(1-p))^-shape-1)/shape)
  return(level)
}
return_level_20 = lapply(monthly_fits, return_level) # 20 for testing
return_level_100 = lapply(monthly_fits, return_level, period=100)
return_level_50 = lapply(monthly_fits, return_level, period=50)
plot(unlist(return_level_100),main="100 Year Return level, estimated with point process", xlab="Month",ylab="Return Level")

plot(unlist(return_level_50),main="100 Year Return level, estimated with point process", xlab="Month",ylab="Return Level")

TODO: estimat with mcmc Assuming that we have the posterior densities of the markov chains, call theis function to plot return level histograms

return_level_mcmc = function(posterior,period=20){
  u = mc.quant(posterior,p=1-1/period,lh="gev")
  label_mcmc_rl = sprintf("%s Year return level",period)
  hist(u,nclass=20,prob=T,xlab=label_mcmc_rl)
}
LS0tCnRpdGxlOiAiVGh1bmRlcnN0cm9tIEFuYWx5c2lzIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCioqSW1wb3J0cyoqCmBgYHtyfQpsaWJyYXJ5KGV2ZCk7IApsaWJyYXJ5KGV2ZGJheWVzKTsKbGlicmFyeShjb2RhKTsKbGlicmFyeShpc21ldik7CmxpYnJhcnkoeHRzKTsKCmBgYAoKCioqTG9hZGluZyB0aGUgRGF0YSoqCgpgYGB7cn0KbG9hZCgiLi4vZGF0YS9DQVBFX01pbmRlcl9SeWNoZW5lcl9NYWxzb3QuUkRhdGEiKQpsb2FkKCIuLi9kYXRhL05JTk8zNC5SRGF0YSIpCmxvYWQoIi4uL2RhdGEvU1JIX01pbmRlcl9SeWNoZW5lcl9NYWxzb3QuUkRhdGEiKQpgYGAKCioqR2VuZXJhdGUgUFJPRCoqCmBgYHtyfQpwcm9kID0gc3FydChjYXBlKSpzcmgKYGBgCgoKCgoqKiBDcmVhdGUgVGltZSBTZXJpZXMgT2JqZWN0cyAqKgpgYGB7cn0KZGF0ZXMgPC0gc2VxLkRhdGUoYXMuRGF0ZSgiMTk3OS0xLTEiKSwgYXMuRGF0ZSgiMjAxNS0xMi0zMSIpLCBieT0xKQpmZWIyOWl4IDwtIGZvcm1hdChhcy5EYXRlKGRhdGVzKSwgIiVtLSVkIikgPT0gIjAyLTI5IgpkYXRlcyA8LSBkYXRlc1shZmViMjlpeF0KCnByb2RfdHMgPSB4dHMocHJvZCwgb3JkZXIuYnkgPSByZXAoZGF0ZXMsIGVhY2g9OCkpCmBgYAoKCgoqKkJlZ2lubmluZyB0aGUgQW5hbHlzaXMqKgpgYGB7cn0KbW9udGhfbmFtZXMgPSBjKCJqYW4iLCJmZWIiLCJtYXIiLCJhcHIiLCJtYXkiLCJqdW4iLCJqdWwiLCJhdWciLCJzZXAiLCJvY3QiLCJub3YiLCJkZWMiKQpnZXRfbW9udGhseSA9IGZ1bmN0aW9uKHgpIHsKICBvdXRwdXQgPSBsaXN0KCkKICBsZW4gPSBucm93KHgpCiAgCiAgIyBHZXQgbW9udGggb2YgZWxlbWVudAogIG1vbnRoID0gdGltZSh4KQogIG1vbnRoID0gZ3N1YigiLi4uLi0iLCAiIiwgbW9udGgpCiAgbW9udGggPSBnc3ViKCItLi4iLCAiIiwgbW9udGgpCiAgbW9udGhsaXN0ID0gdW5pcXVlKG1vbnRoKQogIGZvciAoaSBpbiAxOjEyKSB7CiAgICBvdXRwdXRbW2ldXSA9IHhbbW9udGggPT0gbW9udGhsaXN0W2ldLF0KICB9CiAgbmFtZXMob3V0cHV0KSA9IG1vbnRoX25hbWVzCiAgcmV0dXJuKG91dHB1dCkKfQoKbW9udGhseV9tYXggPSBnZXRfbW9udGhseShhcHBseS5tb250aGx5KHByb2RfdHMsIG1heCkpCnIgPSAyCnJfbW9udGhseV9tYXggPSBnZXRfbW9udGhseShhcHBseS5tb250aGx5KHByb2RfdHMsIGZ1bmN0aW9uKHgpIG9yZGVyKHgsIGRlY3JlYXNpbmc9VClbMTpyXSkpCmBgYAoKCmBgYHtyfQojIGdldCB0aGUgbW9udGhseSBtYXhpbWEgb2YgcHJvZAptMSA9IGFzLmRhdGEuZnJhbWUoYXBwbHkubW9udGhseShwcm9kX3RzLCBtYXgpKSRWMTsKIyBwcm9kdWNlIHRoZSBndW1iZWwgcXEgcGxvdApndW1iZWxfcXEgPSBmdW5jdGlvbih4LCB0aXRsZT0iIikgcXFwbG90KHFndW1iZWwoYygxOmxlbmd0aCh4KSkvKGxlbmd0aCh4KSsxKSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtYWluID0gcGFzdGUoIkd1bWJlbGwgUS1RIFBsb3QiLCB0aXRsZSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeGxhYiA9ICJUaGVvcmV0aWNhbCBRdWFudGlsZXMiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5bGFiID0gIlNhbXBsZSBRdWFudGlsZXMiKSA7CgpndW1iZWxfcXEobTEpCgojcXFwbG90KHFndW1iZWwoYygxOmxlbmd0aChtMSkpLyhsZW5ndGgobTEpKzEpKSxtMSxtYWluID0gIkd1bWJlbGwgUS1RIFBsb3QiLHhsYWIgPSAiVGhlb3JldGljYWwgUXVhbnRpbGVzIiwgeWxhYiA9ICJTYW1wbGUgUXVhbnRpbGVzIikgOwpgYGAKVGhlIHFxIHBsb3QgZG9lc24ndCBmaXQgdmVyeSB3ZWxsLCBlc3BlY2lhbGx5IGluIHRoZSBsb3dlciB0YWlsLiBUaGlzIGlzIGxpa2VseSBkdWUKdG8gc2Vhc29uYWwgZGVwZW5kZW5jZS4KCioqRml0dGluZyBHRVYgdG8gdGhlIGVudGlyZSBkYXRhKioKYGBge3J9CiMgZml0IGdldmQgd2l0aCBNTEUgYW5kIHByb2R1Y2UgZGlhZ25vc3RpYyBwbG90cwpmaXRtYXguTUxFPC1mZ2V2KG0xLG1ldGhvZD0iTmVsZGVyLU1lYWQiKQpwYXIobWZyb3c9YygyLDIpKQpmaXRtYXguTUxFCnBsb3QoZml0bWF4Lk1MRSkKYGBgClBvb3IgZml0LCBwcm9iYWJseSBiZWNhdXNlIHRoZSBkaXN0cmlidXRpb24gaXNuJ3Qgc3RhdGlvbmFyeS4gVGhpcyBpcyBlc3BlY2lhbGx5IAp2aXNpYmxlIGluIHRoZSBQcm9iYWJpbGl0eSBwbG90LCBpbiB3aGljaCB0aGUgY29uZmlkZW5jZSBiYW5kIGlzIHN1cnBhc3NlZCwgCmluZGljYXRpbmcgYSBwb29yIGZpdC4KCgpgYGB7cn0KIyBmaXQgZ2V2ZCB3aXRoIEJheWVzaWFuIFRlY2huaXF1ZXMKIyB1c2UgdGhlIHByZXZpb3VzIG91dHB1dHMgKHJvdW5kZWQpIGFzIGluaXRpYWwgdmFsdWVzICh1c2UgYSBkaWZmZXJlbnQgc2hhcGUpCmluaXQ8LWMoMS42ZTQsNGUzLDAuMSkKIyBhcmJpdHJhcnkgcHJpb3JzCm1hdCA8LSBkaWFnKGMoMTAwMDAsMTAwMDAsMTAwKSkgCnBuIDwtIHByaW9yLm5vcm0obWVhbj1jKDAsMCwwKSxjb3Y9bWF0KQojIHByb3Bvc2FsIHN0YW5kYXJkIGRldmlhdGlvbiAoZm91bmQgYnkgdHJpYWwgYW5kIGVycm9yIHRvIGdldCAyMCU8YWNjZXB0YW5jZSByYXRlPDQwJSkKcHNkPC1jKDUwMCwwLjAzLDAuMDIpCiMgZHJhdyAzayBzYW1wbGVzIGZyb20gbWFya292IGNoYWluCm5pdCA9IDMwMDAwCnRoaW5uaW5nID0gMzAwCnBvc3Q8LXBvc3RlcmlvcihuaXQsIGluaXQ9aW5pdCxwcmlvcj1wbixsaD0iZ2V2IixkYXRhPW0xLHBzZD1wc2QpCiMgZGlhZ25vc3RpYyBwbG90cwpNQ01DPC1tY21jKHBvc3RbMTpuaXQgJSUgdGhpbm5pbmcgPT0gMCwgXSwgdGhpbj0zMDApIApwbG90KE1DTUMpIAphdHRyKG1jbWMocG9zdCksImFyIikKCmBgYAoKCmBgYHtyfQojTUNNQ19zdGFiIDwtIG1jbWMocG9zdCwgdGhpbj01MCwgc3RhcnQ9MTAwMCkKYWNmKG1jbWMocG9zdFsxOm5pdCAlJSB0aGlubmluZyA9PSAwLCBdKSkKYGBgCldlIG9ic2VydmUgdGhhdCB0aGVyZSBzZWVtIHRvIGJlIG5vIHN1YnN0YW50aWFsIGlzc3VlcyBmcm9tIHRoZSBhdXRvY29ycmVsYXRpb24gCnBsb3RzLiAKCmBgYHtyfQphcHBseShtY21jKHBvc3RbMTpuaXQgJSUgdGhpbm5pbmcgPT0gMCwgXSksMixtZWFuKQphcHBseShtY21jKHBvc3RbMTpuaXQgJSUgdGhpbm5pbmcgPT0gMCwgXSksMixzZCkKCmBgYAoKKiogRml0IHdpdGggTUxFIGZvciBtb250aHMgc2VwYXJhdGVseSoqCmBgYHtyfQojbW9udGhseV9maXRzID0gbGFwcGx5KG1vbnRobHlfbWF4LCAKIyAgICAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSBmZ2V2KGRhdGEuZnJhbWUoeClbLDFdLCBtZXRob2Q9Ik5lbGRlci1NZWFkIikpCm1vbnRobHlfZml0cyA9IGxpc3QoKQplcnJvcl9jYXNlcyA9IGMoOSwgMTIpCmZvciAoaSBpbiAxOmxlbmd0aChtb250aGx5X21heCkpIHsKICBwcmludChpKQogIAogIGlmIChpICVpbiUgZXJyb3JfY2FzZXMpIHsKICAgIG1vbnRobHlfZml0c1tbaV1dID0gZmdldihhcy5kYXRhLmZyYW1lKG1vbnRobHlfbWF4W1tpXV0pJFYxLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtZXRob2QgPSAiTmVsZGVyLU1lYWQiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0ZC5lcnIgPSBGQUxTRSkKICB9CiAgZWxzZSB7CiAgICBtb250aGx5X2ZpdHNbW2ldXSA9IGZnZXYoYXMuZGF0YS5mcmFtZShtb250aGx5X21heFtbaV1dKSRWMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gIk5lbGRlci1NZWFkIikKICB9Cn0KbmFtZXMobW9udGhseV9maXRzKSA9IG5hbWVzKG1vbnRobHlfbWF4KQoKYGBgCgpgYGB7cn0KZ3VtYmVsX3FxKGRhdGEuZnJhbWUobW9udGhseV9tYXhbWzldXSlbLDFdLCAiU2VwdGVtYmVyIikKZ3VtYmVsX3FxKGRhdGEuZnJhbWUobW9udGhseV9tYXhbWzEyXV0pWywxXSwgIkRlY2VtYmVyIikKYGBgCmBgYHtyfQpnZXRfc2UgPSBmdW5jdGlvbih4LCBpeCkgewogIGlmIChpcy5udWxsKHgkc3RkLmVycikpIDAKICBlbHNlIHgkc3RkLmVycltpeF0KfQptbGVfbG9jID0gdW5saXN0KGxhcHBseShtb250aGx5X2ZpdHMsIGZ1bmN0aW9uKHgpIHgkZXN0aW1hdGVbMV0pKQptbGVfbG9jX3NlID0gdW5saXN0KGxhcHBseShtb250aGx5X2ZpdHMsIGdldF9zZSwgMSkpCm1sZV9zY2FsZSA9IHVubGlzdChsYXBwbHkobW9udGhseV9maXRzLCBmdW5jdGlvbih4KSB4JGVzdGltYXRlWzJdKSkKbWxlX3NjYWxlX3NlID0gdW5saXN0KGxhcHBseShtb250aGx5X2ZpdHMsIGdldF9zZSwgMikpCm1sZV9zaGFwZSA9IHVubGlzdChsYXBwbHkobW9udGhseV9maXRzLCBmdW5jdGlvbih4KSB4JGVzdGltYXRlWzNdKSkKbWxlX3NoYXBlX3NlID0gdW5saXN0KGxhcHBseShtb250aGx5X2ZpdHMsIGdldF9zZSwgMykpCmBgYAoKYGBge3J9CnBsb3Rfd19lcnIgPSBmdW5jdGlvbih4LCB5LCBzZSwgdGl0bGUgPSBOVUxMKSB7CiAgbWF4X2l4ID0gd2hpY2gubWF4KHkpCiAgbWluX2l4ID0gd2hpY2gubWluKHkpCiAgcGxvdCh4LCB5LAogICAgICAgeWxpbSA9IGMoeVttaW5faXhdIC0gc2VbbWluX2l4XSwgeVttYXhfaXhdICsgc2VbbWF4X2l4XSksCiAgICAgICBtYWluID0gdGl0bGUpCiAgYXJyb3dzKHgseS1zZSx4LHkrc2UsIGNvZGU9MywgbGVuZ3RoPTAuMDIsIGFuZ2xlID0gOTApCn0KcGxvdF93X2VycigxOjEyLCBtbGVfbG9jLCBtbGVfbG9jX3NlLCAiTG9jYXRpb24gUGFyYW1ldGVyIGFzIEVzdGltYXRlZCB3aXRoIExpa2VsaWhvb2QiKQpwbG90X3dfZXJyKDE6MTIsIG1sZV9zY2FsZSwgbWxlX3NjYWxlX3NlLCAiU2NhbGUgUGFyYW1ldGVyIGFzIEVzdGltYXRlZCB3aXRoIExpa2VsaWhvb2QiKQpwbG90X3dfZXJyKDE6MTIsIG1sZV9zaGFwZSwgbWxlX3NoYXBlX3NlLCAiU2hhcGUgUGFyYW1ldGVyIGFzIEVzdGltYXRlZCB3aXRoIExpa2VsaWhvb2QiKQoKYGBgCgoqKiBGaXQgd2l0aCBCYXllc2lhbiBNZXRob2RzIGZvciBtb250aHMgc2VwYXJhdGVseSoqCmBgYHtyfQoKCiMgRml0cyBHRVYgZGlzdHJpYnV0aW9uIHdpdGggYmF5ZXNpYW4gbWV0aG9kIGZvciBnaXZlbiBwYXJhbWV0ZXJzCmJheWVzX2ZpdHRlciA9IGZ1bmN0aW9uKHgsIAogICAgICAgICAgICAgICAgICAgICAgICBpbml0ID0gYygxZTMsIDFlMywgMC4xKSwgIyBJbml0aWFsIHZhbHVlcwogICAgICAgICAgICAgICAgICAgICAgICBtYXQgPSBkaWFnKGMoMTAwMDAsMTAwMDAsMTAwKSksCiAgICAgICAgICAgICAgICAgICAgICAgIHBzZCA9IGMoNTAwLDAuMSwwLjEpLCAjIFByb3Bvc2VkIFNEZXYKICAgICAgICAgICAgICAgICAgICAgICAgbml0ID0gMzAwMCwgIyBOYiBJdGVyYXRpb25zCiAgICAgICAgICAgICAgICAgICAgICAgIHRoaW5uaW5nID0gNTAsICMgVGhpbm5pbmcgRmFjdG9yCiAgICAgICAgICAgICAgICAgICAgICAgIGRvX2RpYWduID0gRkFMU0UsICMgQm9vbCB3aGV0aGVyIHRvIHNob3cgZGlhZ25vc3RpYyBwbG90cwogICAgICAgICAgICAgICAgICAgICAgICBkb19hdXRvcmVnID0gRkFMU0UsICMgQm9vbCB3aGV0aGVyIHRvIHNob3cgYXV0b3JlZyBwbG90cwogICAgICAgICAgICAgICAgICAgICAgICBzZWVkID0gNDIgIyBTZWVkIAogICAgICAgICAgICAgICAgICAgICAgICApIHsKICBzZXQuc2VlZChzZWVkKSAgICAgICAgICAgICAgICAKICBwbiA9IHByaW9yLm5vcm0obWVhbj1jKDAsMCwwKSxjb3Y9bWF0KQogIHBvc3QgPSBwb3N0ZXJpb3Iobml0LCBpbml0PWluaXQsIHByaW9yPXBuLCBsaD0iZ2V2IiwgZGF0YT14LCBwc2Q9cHNkKQogIAogIGlmKGRvX2RpYWduKSB7CiAgICBNQ01DPC1tY21jKHBvc3RbMTpuaXQgJSUgdGhpbm5pbmcgPT0gMCwgXSkgCiAgICBwbG90KE1DTUMpIAogIH0KICBpZihkb19hdXRvcmVnKSB7CiAgICBhY2YobWNtYyhwb3N0WzE6bml0ICUlIHRoaW5uaW5nID09IDAsIF0pKQogIH0KICBsaXN0KHBvc3RlcmlvciA9IHBvc3QsIAogICAgICAgYWNjZXB0YW5jZV9yYXRlID0gYXR0cihtY21jKHBvc3QpLCJhciIpKQp9CgoKCiMgSXRlcmF0aXZlbHkgZml0cyBHRVYgd2l0aCBiYXllc2lhbiBtZXRob2RzLCB1bnRpbCB0aGUgZml0IGhhcyAKIyBhY2NlcHRhYmxlIGFjY2VwdGFuY2UgcmF0ZXMgKGkuZS4gMC4yIDwgQVIgPCAwLjQpLiBJZiB0aGUgQVIgaXMgdG9vIGhpZ2gsIAojIHRoZSBwcm9wb3NlZCBTRGV2IGZvciB0aGUgcGFyYW1ldGVyIGlzIG11bHRpcGxpZWQgd2l0aCAxLjUuIElmIGl0J3MgdG9vIHNtYWxsLAojIHRoZSBwcm9wb3NlZCBTRGV2IGlzIGRpdmlkZWQgYnkgMi4gVGhpcyBpcyByZXBlYXRlZCB1bnRpbCB0aGUgYWNjZXB0YW5jZSByYXRlCiMgaXMgZ29vZCBmb3IgYWxsIHBhcmFtZXRlcnMsIG9yIG1heF9pdCBpcyByZWFjaGVkLiBUaGVuLCBhIGZpbmFsIG1vZGVsIGlzIGZpdHRlZAojIHdpdGggbW9yZSBpdGVyYXRpb25zLiAKYmF5ZXNfZml0dGVyX3BhcmFtX3NlYXJjaCA9IGZ1bmN0aW9uKHgsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBwc2RfaW5pdCA9IGMoNTAwLDAuMSwwLjEpLCAjIEluaXRpYWwgcHJvcG9zZWQgU0RldgogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbml0X2Z1bGwgPSAzMDAwLCAjIE5iIEl0ZXJhdGlvbnMgZm9yIGZpbmFsIG1vZGVsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBuaXRfc2VhcmNoID0gMTUwLCAjIE5iIEl0ZXJhdGlvbnMgZm9yIHBhcmFtIHNlYXJjaAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZG9fZGlhZ24gPSBGQUxTRSwgIyBCb29sIHdoZXRoZXIgdG8gc2hvdyBkaWFnbm9zdGljIHBsb3RzCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBkb19hdXRvcmVnID0gRkFMU0UsICMgQm9vbCB3aGV0aGVyIHRvIHNob3cgYXV0b3JlZyBwbG90cwogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWF4X2l0ID0gMjAsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAuLi4gIyBBZGRpdGlvbmFsIHBhcmFtcyBwYXNzZWQgdG8gYmF5ZXNfZml0dGVyCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICApIAp7CiAgIyBJdGVyYXRlIHVudGlsIGRlc2lyZWQgYWNjZXB0YW5jZSByYXRlIGlzIG9idGFpbmVkCiAgY29udCA9IFRSVUUKICBwc2QgPSBwc2RfaW5pdAogIGl0ID0gMAogIHdoaWxlKGNvbnQpIHsKICAgIGl0ID0gaXQrMQogICAgaWYgKGl0ID4gbWF4X2l0KSB7CiAgICAgIHdhcm5pbmcoIlRoZSAiKQogICAgfQogICAgZml0ID0gYmF5ZXNfZml0dGVyKHgsIHBzZD1wc2QsIG5pdD1uaXRfc2VhcmNoLCBkb19kaWFnbj1GQUxTRSwgCiAgICAgICAgICAgICAgICAgICAgICAgZG9fYXV0b3JlZz1GQUxTRSwuLi4pCiAgICBhY2NfcmF0ZXMgPSBmaXQkYWNjZXB0YW5jZV9yYXRlWzEsIDE6M10KICAgIAogICAgdG9vX2hpZ2ggPSBhY2NfcmF0ZXMgPiAuNAogICAgdG9vX2xvdyA9IGFjY19yYXRlcyA8IC4yCiAgICAKICAgIGlmIChhbGwoIXRvb19oaWdoKSAmJiBhbGwoIXRvb19sb3cpKSB7ICMgQWxsIGFjY2VwdGFuY2UgcmF0ZXMgbGllIHdpdGhpbiB0aHJlc2hvbGQKICAgICAgY29udD1GQUxTRQogICAgfSBlbHNlIGlmIChpdCA+IG1heF9pdCkgeyAjIG1heF9pdCBpcyByZWFjaGVkCiAgICAgIGNvbnQ9RkFMU0UKICAgICAgd2FybmluZygibWF4X2l0IHdhcyByZWFjaGVkIikKICAgIH0gZWxzZSB7ICMgQ29ycmVjdCB2YWx1ZXMgd2hpY2ggaGF2ZSB3cm9uZyB0aHJlc2hvbGQKICAgICAgcHNkW3Rvb19oaWdoXSA9IHBzZFt0b29faGlnaF0gKiAxLjUKICAgICAgcHNkW3Rvb19sb3ddID0gcHNkW3Rvb19sb3ddIC8gMgogICAgfQogIH0KICAKICAjIEZpdCBmaW5hbCBtb2RlbAogIGJheWVzX2ZpdHRlcih4LCBwc2Q9cHNkLCBuaXQ9bml0X2Z1bGwsIGRvX2RpYWduPWRvX2RpYWduLCAKICAgICAgICAgICAgICAgZG9fYXV0b3JlZz1kb19hdXRvcmVnLCAuLi4pCiAgCn0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIAoKbW9udGhseV9iYXllc19maXQgPSBsYXBwbHkobW9udGhseV9tYXgsIGJheWVzX2ZpdHRlcl9wYXJhbV9zZWFyY2gsIGRvX2RpYWduID0gVFJVRSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgIGRvX2F1dG9yZWcgPSBGQUxTRSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgcHNkID0gYyg1MDAsMC4zLDAuMyksIG5pdF9mdWxsPTMwMDAwLCBuaXRfc2VhcmNoID0gMzAwMDAsCiAgICAgICAgICAgICAgICAgICAgICAgICAgIHRoaW5uaW5nID0gMzAwKQphY2NlcHRhbmNlX3JhdGVzID0gbGFwcGx5KG1vbnRobHlfYmF5ZXNfZml0LCBmdW5jdGlvbih4KSB4JGFjY2VwdGFuY2VfcmF0ZVsxLF0pCnByaW50KGFjY2VwdGFuY2VfcmF0ZXMpCmJheWVzX3BhcmFtcyA9IGxhcHBseShtb250aGx5X2JheWVzX2ZpdCwgZnVuY3Rpb24oeCkgYXBwbHkoeCRwb3N0ZXJpb3IsIDIsIG1lYW4pKQpiYXllc19zdGRlcnIgPSBsYXBwbHkobW9udGhseV9iYXllc19maXQsIGZ1bmN0aW9uKHgpIGFwcGx5KHgkcG9zdGVyaW9yLCAyLCBzZGV2KSkKCgpgYGAKCgpUT0RPLT4gUiBsYXJnZXN0IGZpdAoKCgoqKlBBUlQgMioqCkZpcnN0LCB3ZSBjaGVjayBpZiB0aGUgbG9jYXRpb24gcGFyYW1ldGVyIGRlcGVuZHMgb24gdGltZSB1c2luZyBhIGxpa2VseWhvb2QgcmF0aW8gdGVzdApgYGB7cn0KI21vbnRobHlfZml0cyA9IGxhcHBseShtb250aGx5X21heCwgCiMgICAgICAgICAgICAgICAgICAgICAgZnVuY3Rpb24oeCkgZmdldihkYXRhLmZyYW1lKHgpWywxXSwgbWV0aG9kPSJOZWxkZXItTWVhZCIpKQpyYXRpb3MgPSBsaXN0KCkKdHJlbmQgPSAxOmxlbmd0aChhcy5kYXRhLmZyYW1lKG1vbnRobHlfbWF4W1sxMl1dKSRWMSkKdHJlbmQgPSAodHJlbmQtbWVhbih0cmVuZCkpL3NkKHRyZW5kKSAjIHNjYWxlIGFuZCBjZW50ZXIgY292YXJpYXRlcyBhcyByZWNvbW1lbmRlZAplcnJvcl9jYXNlcyA9IGMoOSwgMTIpCmZvciAoaSBpbiAxOmxlbmd0aChtb250aGx5X21heCkpIHsKICBwcmludChpKQogIAogIGZpdF9jb25zdCA9IGZnZXYoYXMuZGF0YS5mcmFtZShtb250aGx5X21heFtbaV1dKSRWMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gIk5lbGRlci1NZWFkIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzdGQuZXJyID0gRkFMU0UpCiAgZml0X2RlcGVuZGFudCA9IGZnZXYoYXMuZGF0YS5mcmFtZShtb250aGx5X21heFtbaV1dKSRWMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gIk5lbGRlci1NZWFkIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBuc2xvYyA9IHRyZW5kLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0ZC5lcnIgPSBGQUxTRSkKICAKICByYXRpb3NbW2ldXSA9IGZpdF9jb25zdCRkZXYtZml0X2RlcGVuZGFudCRkZXYgCn0KbmFtZXMocmF0aW9zKSA9IG5hbWVzKG1vbnRobHlfbWF4KQpjaGlfOTVsZXZlbCA9IHFjaGlzcSgxLTAuMDUvMTIsMSkKCnBsb3QodW5saXN0KHJhdGlvcyksbWFpbj0iOTUlIGNvbmZpZGVuY2UgdGVzdCBmb3IgdGltZSBpbmRlcGVuZGFuY2UsIEJvbmZlcnJvbmkgbXVsdGlwbGUgVGVzdGluZyIsIHhsYWI9Ik1vbnRoIix5bGFiPSJMaWtlbHlob29kIFJhdGlvIFN0YXRpc3RpYyIpCmFibGluZShhPWNoaV85NWxldmVsLGI9MCxjb2w9InJlZCIpCgpgYGAKCk5vdywgbGV0J3MgY2hlY2sgZm9yIGluZGVwZW5kYW5jZSBmcm9tIEVOU08KYGBge3J9CiNtb250aGx5X2ZpdHMgPSBsYXBwbHkobW9udGhseV9tYXgsIAojICAgICAgICAgICAgICAgICAgICAgIGZ1bmN0aW9uKHgpIGZnZXYoZGF0YS5mcmFtZSh4KVssMV0sIG1ldGhvZD0iTmVsZGVyLU1lYWQiKSkKcmF0aW9zID0gbGlzdCgpCiMgc3BsaXQgbmlubyBkYXRhIGludG8gbW9udGhzCm4gPSBuaW5vMzQKZGltKG4pPWMoMTIsbGVuZ3RoKGFzLmRhdGEuZnJhbWUobW9udGhseV9tYXhbWzEyXV0pJFYxKSkKZXJyb3JfY2FzZXMgPSBjKDksIDEyKQpmb3IgKGkgaW4gMTpsZW5ndGgobW9udGhseV9tYXgpKSB7CiAgcHJpbnQoaSkKICB0cmVuZCA9IG5baSxdCiAgdHJlbmQgPSAodHJlbmQtbWVhbih0cmVuZCkpL3NkKHRyZW5kKSAjIHNjYWxlIGFuZCBjZW50ZXIgY292YXJpYXRlcyBhcyByZWNvbW1lbmRlZAogIGZpdF9jb25zdCA9IGZnZXYoYXMuZGF0YS5mcmFtZShtb250aGx5X21heFtbaV1dKSRWMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gIk5lbGRlci1NZWFkIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzdGQuZXJyID0gRkFMU0UpCiAgZml0X2RlcGVuZGFudCA9IGZnZXYoYXMuZGF0YS5mcmFtZShtb250aGx5X21heFtbaV1dKSRWMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gIk5lbGRlci1NZWFkIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBuc2xvYyA9IHRyZW5kLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0ZC5lcnIgPSBGQUxTRSkKICAKICByYXRpb3NbW2ldXSA9IGZpdF9jb25zdCRkZXYtZml0X2RlcGVuZGFudCRkZXYgCn0KbmFtZXMocmF0aW9zKSA9IG5hbWVzKG1vbnRobHlfbWF4KQpjaGlfOTVsZXZlbCA9IHFjaGlzcSgxLTAuMDUvMTIsMSkKCnBsb3QodW5saXN0KHJhdGlvcyksbWFpbj0iOTUlIGNvbmZpZGVuY2UgdGVzdCBmb3IgaW5kZXBlbmRhbmNlIGZyb20gRU5TTywgQm9uZmVycm9uaSBtdWx0aXBsZSB0ZXN0aW5nIiwgeGxhYj0iTW9udGgiLHlsYWI9Ikxpa2VseWhvb2QgUmF0aW8gU3RhdGlzdGljIikKYWJsaW5lKGE9Y2hpXzk1bGV2ZWwsYj0wLGNvbD0icmVkIikKCmBgYAoKQW5vdGhlciBtZXRob2QgaXMgdGhlIGNoaSBwbG90OgpgYGB7cn0KIyBwbG90IHRoZSBjaGkgcGxvdCBmb3IgZGVwZW5kYW5jZSBhbmFseXNpcwp0cmVuZCA9IDE6bGVuZ3RoKGFzLmRhdGEuZnJhbWUobW9udGhseV9tYXhbWzEyXV0pJFYxKQpuID0gbmlubzM0CmRpbShuKT1jKDEyLGxlbmd0aChhcy5kYXRhLmZyYW1lKG1vbnRobHlfbWF4W1sxMl1dKSRWMSkpCgpmb3IgKGkgaW4gMTpsZW5ndGgobW9udGhseV9tYXgpKSB7CiAgbmlubyA9IG5baSxdCiAgbV9kYXRhPWFzLmRhdGEuZnJhbWUobW9udGhseV9tYXhbW2ldXSkkVjEKICBkYXQubTFfbW9udGggPSBjYmluZChtX2RhdGEsdHJlbmQpOwogIGRhdC5tMV9uaW5vID0gY2JpbmQobV9kYXRhLG5pbm8pOwogIHBhcihtZnJvdz1jKDIsMikpCiAgY2hpcGxvdChkYXQubTFfbW9udGgsbWFpbjEgPSAiQ2hpIFBsb3QgVGltZSIsbWFpbjIgPSAiQ2hpIEJhciBQbG90IFRpbWUiKTsKICBjaGlwbG90KGRhdC5tMV9uaW5vLG1haW4xID0gIkNoaSBQbG90IEVOU08iLG1haW4yID0gIkNoaSBCYXIgUGxvdCBFTlNPIik7Cn0KYGBgCgoKCioqUEFSVCAzKioKV2Ugd2lsbCBub3cgYW5hbHlzZSB0ZW1wb3JhbCBjbHVzdGVyaW5nIG9mIGV4dHJlbWVzLiBGb3IgdGhpcywgd2Ugd2lsbCB1c2UgdGhlIGV4aXBsb3QgZnVuY3Rpb24gZnJvbSB0aGUgZXZkIGxpYnJhcnkuCgpgYGB7cn0KIyBkZWZpbmUgYSBmdW5jdGlvbiBmb3IgZ2V0dGluZyB0aGUgZXh0cmVtYWwgaW5kaWNlcyBmb3IgZWFjaCBtb250aCBmb3IgYSBnaXZlbiB0aHJlc2hvbGQKbW9udGhseV9laW5kZXggPC0gZnVuY3Rpb24oZGF0YSwgdGhyZXNob2xkX3AsIHI9MCl7CiAgZWkgPSBsaXN0KCkKICBmb3IgKGkgaW4gMTpsZW5ndGgoZGF0YSkpIHsKICAgIHRocmVzaG9sZCA9IHF1YW50aWxlKGFzLmRhdGEuZnJhbWUoZGF0YVtbaV1dKSRWMSwgdGhyZXNob2xkX3ApCiAgICBlaVtbaV1dPWV4aShhcy5kYXRhLmZyYW1lKGRhdGFbW2ldXSkkVjEsIHRocmVzaG9sZCwgcikKICB9CiAgbmFtZXMoZWkpID0gbmFtZXMoZGF0YSkKCiAgcmV0dXJuKGVpKQp9CgplaSA9IG1vbnRobHlfZWluZGV4KGdldF9tb250aGx5KHByb2RfdHMpLCAwLjk1KQpwbG90KHVubGlzdChlaSksIG1haW49IkV4dHJlbWFsIEluZGV4IGJ5IE1vbnRoLCA5NSUtUXVhbnRpbGUgYXMgVGhyZXNob2xkIiwgeGxhYj0iTW9udGgiLCB5bGFiPSJFeHRyZW1hbCBpbmRleCIpCmBgYAoKV2Ugb2JzZXJ2ZSB0aGF0IHRoZSBleHRyZW1hbCBpbmRleCBpcyB+MC4yNS1+MC40NSwgd2UgY2FuIHRoZXJlZm9yZSBjb25jbHVkZSB0aGF0IHdlIGhhdmUgc3Ryb25nIGRlcGVuZGFuY2Ugb2YgZXh0cmVtZXMsIHdpdGggdGhlIGxpbWl0aW5nIG1lYW4gY2x1c3RlciBzaXplIGJlaW5nIHJvdWdobHkgZnJvbSAyIHRvIDQuIFRoZSBjbHVzdGVyaW5nIGhhcyBubyBlZmZlY3QgZm9yIGVzdGltYXRlcnMgYmFzZWQgb24gdGhlIChtb250aGx5KSBtYXhpbXVtLCBidXQgdGhlIHIgbGFyZ2VzdCBlc3RpbWF0b3IgaXMgaW5mbHVlbmNlZCBieSBpdC4KCgoqKlBBUlQgNCoqCkZpcnN0LCBsZXQncyBlc3RpbWF0ZSB0aGUgMTAgeWVhciByZXR1cm4gbGV2ZWwgdXNpbmcgcG9pbnQgcHJvY2VzcyBhcHByb2FjaApgYGB7cn0KbW9udGhseV9maXRzX3BwID0gbGlzdCgpCm1vbnRobHlfZGF0YSA9IGdldF9tb250aGx5KHByb2RfdHMpCiNlcnJvcl9jYXNlcyA9IGMoMSwyLDMsNCw1LDYsNyw4LDksMTAsMTEsMTIpCm1vbnRoX2RheXMgPSBjKDMxLDI4LDMxLDMwLDMxLDMwLDMxLDMxLDMwLDMxLDMwLDMxKQpmb3IgKGkgaW4gMTpsZW5ndGgobW9udGhseV9tYXgpKSB7CiAgcHJpbnQoaSkKICB0aHJlc2hvbGQgPSBxdWFudGlsZShhcy5kYXRhLmZyYW1lKG1vbnRobHlfZGF0YVtbaV1dKSRWMSwgMC45NSkKICAKICBpZiAoaSAlaW4lIGVycm9yX2Nhc2VzKSB7CiAgICBtb250aGx5X2ZpdHNfcHBbW2ldXSA9IGZwb3QoYXMuZGF0YS5mcmFtZShtb250aGx5X2RhdGFbW2ldXSkkVjEsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgdGhyZXNob2xkID0gdGhyZXNob2xkLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIG1vZGVsPSJwcCIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbnBwID0gbW9udGhfZGF5c1tpXSo4LAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIGNtYXggPSBUUlVFLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHIgPSAxLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0ZC5lcnIgPSBGQUxTRSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtZXRob2QgPSAiTmVsZGVyLU1lYWQiKQogIH0KICBlbHNlIHsKICAgIG1vbnRobHlfZml0c19wcFtbaV1dID0gZnBvdChhcy5kYXRhLmZyYW1lKG1vbnRobHlfZGF0YVtbaV1dKSRWMSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICB0aHJlc2hvbGQgPSB0aHJlc2hvbGQsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbW9kZWw9InBwIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBucHAgPSBtb250aF9kYXlzW2ldKjgsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgY21heCA9IFRSVUUsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgciA9IDEsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gIk5lbGRlci1NZWFkIikKICB9Cn0KbmFtZXMobW9udGhseV9maXRzX3BwKSA9IG5hbWVzKG1vbnRobHlfZGF0YSkKZm9yKGkgaW4gMToxMil7CiAgcGFyKG1mcm93PWMoMiwyKSkgCiAgcGxvdChtb250aGx5X2ZpdHNfcHBbW2ldXSkKfQoKYGBgClRoZSBmaXQgaW4gZmVicnVhcnkgaGFzIGNvbXBsZXRlbHkgZmFpbGVkLCBhbmQgdGhlIG90aGVycyBhcmUgbm90IHZlcnkgZ29vZCBlaXRoZXIKCldlIHdpbGwgc3RpbGwgZXN0aW1hdGUgdGhlIHJldHVybiBsZXZlbHM6CmBgYHtyfQpyZXR1cm5fbGV2ZWwgPSBmdW5jdGlvbih4LHBlcmlvZD0yMCl7CiAgcCA9IDEvcGVyaW9kCiAgbG9jID0geCRlc3RpbWF0ZVtbMV1dCiAgc2NhbGUgPSB4JGVzdGltYXRlW1syXV0KICBzaGFwZSA9IHgkZXN0aW1hdGVbWzNdXQogIGxldmVsID0gbG9jICsgc2NhbGUqKCgoLWxvZygxLXApKV4tc2hhcGUtMSkvc2hhcGUpCiAgcmV0dXJuKGxldmVsKQp9CnJldHVybl9sZXZlbF8yMCA9IGxhcHBseShtb250aGx5X2ZpdHMsIHJldHVybl9sZXZlbCkgIyAyMCBmb3IgdGVzdGluZwpyZXR1cm5fbGV2ZWxfMTAwID0gbGFwcGx5KG1vbnRobHlfZml0cywgcmV0dXJuX2xldmVsLCBwZXJpb2Q9MTAwKQpyZXR1cm5fbGV2ZWxfNTAgPSBsYXBwbHkobW9udGhseV9maXRzLCByZXR1cm5fbGV2ZWwsIHBlcmlvZD01MCkKcGxvdCh1bmxpc3QocmV0dXJuX2xldmVsXzEwMCksbWFpbj0iMTAwIFllYXIgUmV0dXJuIGxldmVsLCBlc3RpbWF0ZWQgd2l0aCBwb2ludCBwcm9jZXNzIiwgeGxhYj0iTW9udGgiLHlsYWI9IlJldHVybiBMZXZlbCIpCnBsb3QodW5saXN0KHJldHVybl9sZXZlbF81MCksbWFpbj0iNTAgWWVhciBSZXR1cm4gbGV2ZWwsIGVzdGltYXRlZCB3aXRoIHBvaW50IHByb2Nlc3MiLCB4bGFiPSJNb250aCIseWxhYj0iUmV0dXJuIExldmVsIikKCmBgYAoKClRPRE86IGVzdGltYXQgd2l0aCBtY21jCkFzc3VtaW5nIHRoYXQgd2UgaGF2ZSB0aGUgcG9zdGVyaW9yIGRlbnNpdGllcyBvZiB0aGUgbWFya292IGNoYWlucywgY2FsbCB0aGVpcyBmdW5jdGlvbiB0byBwbG90IHJldHVybiBsZXZlbCBoaXN0b2dyYW1zCgpgYGB7cn0KcmV0dXJuX2xldmVsX21jbWMgPSBmdW5jdGlvbihwb3N0ZXJpb3IscGVyaW9kPTIwKXsKICB1ID0gbWMucXVhbnQocG9zdGVyaW9yLHA9MS0xL3BlcmlvZCxsaD0iZ2V2IikKICBsYWJlbF9tY21jX3JsID0gc3ByaW50ZigiJXMgWWVhciByZXR1cm4gbGV2ZWwiLHBlcmlvZCkKICBoaXN0KHUsbmNsYXNzPTIwLHByb2I9VCx4bGFiPWxhYmVsX21jbWNfcmwpCn0KCgpgYGAKCgoKCgoKCgo=